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1.0  SUMMARY 


A  primary  aspect  involved  in  the  navigation  of  a  wide  variety  of  vehicles,  such  as  ground  vehicles 
and  aircraft,  is  the  use  of  inertial  measurement  units.  The  purpose  of  the  inertial  measurement 
unit  is  to  provide  an  internal  measure  of  the  linear  acceleration  and  angular  velocity  of  the  vehicle 
in  order  to  be  used  by  a  navigation  system  in  the  dead-reckoning  of  the  vehicle’s  position  and 
attitude.  Due  to  the  nature  of  sensors,  the  measurement  of  the  linear  acceleration  and  angular 
velocity  are  corrupted  by  several  error  sources,  such  as  bias  and  noise,  amongst  others.  In  order 
to  assess  the  coupled  impact  of  these  error  sources  on  the  navigation  accuracy  for  an  arbitrary 
vehicle,  a  simulation  and  analysis  tool  is  developed  that  can  model  several  inertial  measurement 
unit  mechanizations,  simulate  the  error  sources  encountered  in  the  acquisition  of  measurement 
data,  emulate  the  navigation  software,  and  perform  a  range  of  statistical  analyses  via  Monte  Carlo 
simulation.  The  developed  tool  is  applied  to  a  ground  vehicle  navigation  scenario  in  order  to  assess 
the  performance  of  different  sensors  and  provide  conclusions  as  to  when  each  sensor  should  be  used 
to  achieve  the  best  navigation  performance. 

2.0  INTRODUCTION 

Inertial  measurement  unit  (IMU)  technology  is  used  to  drive  a  dominant  element  of  navigation 
systems  on  vehicles  ranging  from  submarines  to  ground  vehicles,  missiles,  and  spacecraft.  The 
purpose  of  the  IMU  is  to  provide  an  internal  measure  of  the  changes  in  the  linear  and  angular 
motions  of  the  vehicle  due  to  non-gravitational  stimuli  by  using  triads  of  single-axis  accelerometers 
to  measure  the  linear  motion  and  gyros  to  measure  the  angular  motion  via  measurement  of  the 
angular  velocity.  In  doing  so,  the  position,  velocity,  and  attitude  of  the  vehicle  can  be  predicted 
from  a  prescribed  starting  condition  along  with  a  representation  of  the  confidence  in  the  prediction, 
i.e.  an  associated  uncertainty. 

Traditionally,  IMU-based  navigation  is  mechanized  in  one  of  two  ways:  gimbaled/gyro-stabilized 
or  strapdown.  In  the  gimbaled  mechanizations,  the  triad  of  accelerometers  is  placed  on  a  “stable 
table”  at  the  center  of  set  of  interconnected  gimbals  or  at  the  center  of  a  fluidic  stabilization  system. 
In  the  strapdown  mechanization,  the  triad  of  accelerometers  and  gyros  are  fixed  (strapped  down) 
to  the  body  of  the  vehicle,  and  there  is  no  stable  table.  The  different  mechanizations  effectively 
alter  the  way  in  which  the  triad  orientation  is  handled:  for  gimbaled  systems,  mechanical  tracking 
is  used,  whereas  computational  tracking  is  used  for  strapdown  systems.  Common  to  every  system, 
the  outputs  of  the  single-axis  sensors  that  make  up  the  accelerometer  and  gyroscope  are  used  in 
conjunction  with  the  navigation  system  software  to  dead-reckon  (i.e.  using  only  the  IMU)  the 
state  (position,  velocity,  and  attitude)  of  the  vehicle.  Since  the  single-axis  sensors  are  imperfect 
(corrupted  by  measurement  errors),  knowledge  of  the  exact  position,  velocity,  and  attitude  of  the 
vehicle  will  degrade  over  time,  which  leads  to  navigation  errors  that  may  be  described  by  a  statistical 
distribution. 

Typical  analysis  of  IMU-driven  systems  employs  limited  statistical  information  regarding  the 
distribution  of  the  position,  velocity,  and  attitude  of  the  vehicle.  Namely,  the  standard  linear 
covariance  analysis  makes  use  of  only  the  first  two  central  moments  of  the  distribution,  i.e.  the 
mean  and  covariance.  The  covariance  obtained  from  linear  covariance  analysis  is  usually  interpreted 
as  being  representative  of  a  Gaussian  distribution  in  order  to  determine  performance  measures; 
however,  if  the  distribution  is  not  Gaussian,  the  calculation  of  performance  measures  can  lead 
to  erroneous  conclusions.  Additionally,  using  only  the  mean  and  covariance  information  from 
a  Gaussian  distribution  limits  the  type  of  performance  analysis  that  can  be  conducted.  Monte 
Carlo  techniques  present  an  attractive  alternative  that  obviates  erroneous  statistical  conclusions 


1 

Approved  for  public  release;  distribution  is  unlimited. 


obtained  from  a  Gaussian  assumption  resulting  from  linear  covariance  analysis  by  working  with  a 
fully  probabilistic  representation  of  the  navigation  errors. 

The  following  objectives  are  identified  as  necessary  for  the  development  of  a  tool  that  can 
simulate  and  analyze  the  performance  of  IMU-based  navigation  systems: 

1.  Formulate  dynamical  system  representations  of  the  state  evolution  model  for  gimbal-  and 
strap  down-based  systems. 

2.  Develop  a  general  IMU  model  capable  of  simulating  the  corrupted  acceleration  and  angular 
velocity  of  an  IMU  given  true  acceleration  and  angular  velocity  trajectories. 

3.  Identify  relevant  performance  criteria  to  be  assessed,  such  as, 

•  error  probables, 

•  terminal  miss  distance  Monte  Carlo  sample  histograms, 

•  terminal  miss  distance  Monte  Carlo  sample  clouds,  and 

•  position,  velocity,  and  attitude  standard  deviation  over  time. 

4.  Build  a  simulation  environment  to  assess  navigation  performance  which  contains: 

•  a  vehicle  model  (i.e.  gimbal-based  or  strapdown-based  vehicle), 

•  an  IMU  error  model  that  allows  for  user-selectable  inclusion  of  different  error  sources 
and  statistics  of  the  parameterization, 

•  an  IMU-based  dead-reckoning  navigation  model,  and 

•  an  analysis  model  that  allows  for  user-selectable  performance  criteria  to  be  assessed. 

5.  Perform  Monte  Carlo  analyses  and  assess  performance  criteria. 

In  order  to  quantify  the  effects  of  the  IMU  parameters  on  navigation  error,  SAIMUN  (Simulation 
and  Analysis  of  IMU-based  Navigation)  was  developed.  This  tool  utilizes  a  known  truth  trajectory, 
IMU  parameterization  of  interest,  and  initial  conditions  to  calculate  the  performance  criteria.  The 
tool  makes  use  of  Monte  Carlo  analysis  to  propagate  the  position,  velocity,  and  attitude  uncertainty 
forward  in  time.  Each  Monte  Carlo  sample  is  generated  according  to  the  initial  uncertainty  and 
is  then  dead-reckoned  with  corrupted  acceleration  and  angular  velocity  data  according  to  the  IMU 
parameterization.  Statistics  are  collected  from  these  Monte  Carlo  samples  over  time  to  obtain  the 
performance  criteria. 

SAIMUN  was  developed  to  consider  three  different  IMU  mechanizations:  strapdown,  space- 
stabilized,  and  local-level.  The  strapdown  IMU  mechanization  relies  on  microelectromechanical 
systems  (MEMS)  technology  for  the  accelerometer  and  gyroscope  triads.  The  triads  are  stationary 
with  respect  to  the  IMU  case;  thus,  the  attitude  is  dead-reckoned  numerically  using  the  gyroscope 
data  to  provide  the  vehicle  attitude.  The  space-stabilized  and  local-level  IMU  mechanizations  rely 
on  a  mechanical  platform  to  track  a  desired  reference  frame  in  which  the  accelerometer  triad  is 
oriented.  The  local-level  mechanization  tracks  the  North-East-Down  coordinate  frame;  the  space- 
stabilized  mechanization  tracks  an  inertial  coordinate  frame. 

The  remainder  of  this  report  is  organized  as  follows.  Section  3.1  provides  an  overview  of  the 
simulation  environment  that  was  developed.  Section  3.2  details  the  modeling  of  the  error  sources 
which  contribute  to  the  IMU  corruption  and  provides  the  overall  models  for  the  measured  IMU 
outputs  that  are  used  in  the  strapdown,  local-level,  and  space-stabilized  implementations.  In 
Section  3.3,  an  overview  of  the  Allan  variance  analysis  method  and  its  utility  for  analyzing  noise 
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processes  present  in  IMU  data  is  given,  which  is  followed  by  a  detailed  development  of  the  dynamical 
models  used  for  the  strapdown,  local-level,  and  space-stabilized  IMUs  in  Section  3.4.  Section  3.5 
describes  the  Monte  Carlo  analysis  technique  and  its  advantages  over  other  uncertainty  analysis 
techniques.  Section  4.0  presents  the  results  of  a  ground  vehicle  navigation  simulation  in  which  two 
IMUs  are  considered  and  compared  using  the  developed  tool.  Section  5.0  concludes  the  report  with 
a  discussion  as  to  how  the  developed  tool  is  used  to  select  an  appropriate  IMU  for  the  ground 
vehicle  navigation  simulation  as  well  as  a  discussion  on  other  capabilities  that  SAIMUN  offers. 

3.0  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 
3.1  SAIMUN  Overview 

In  order  to  analyze  the  effect  of  inertial  measurement  unit  performance  on  navigation  systems, 
SAIMUN  was  developed.  SAIMUN  utilizes  Monte  Carlo  sampling  to  propagate  the  initial  position, 
velocity,  and  attitude  uncertainty  forward  in  time,  while  accounting  for  the  injection  of  uncertainty 
into  the  position,  velocity,  and  attitude  due  to  error  corruption  present  in  the  IMU.  The  effect  of 
the  following  IMU  error  sources  on  the  position,  velocity,  and  attitude  uncertainty  is  quantified  in 
SAIMUN; 

•  Startup  bias 

•  Bias  instability 

•  Random  walk 

•  Scale  factor  error 

•  Axes  misalignment 

•  Axes  nonorthogonality 

•  Quantization 

The  modeling  of  these  IMU  error  parameters  is  outlined  in  Section  3.2.  In  addition  to  information 
about  the  IMU  error  sources,  the  tool  also  requires  the  following  information: 

•  Tabulated  true  trajectory,  which  is  provided  as  an  input  to  SAIMUN  to  ensure  that  the  tool 
can  be  applied  to  any  vehicle,  ranging  from  ground  vehicles  to  aquatic  vehicles,  air  vehicles, 
and  spacecraft.  The  tabulated  truth  trajectory  consists  of  the  following  information: 

o  inertial  position  vector  of  the  center  of  mass  (CM)  of  the  vehicle, 
o  inertial  velocity  vector  of  the  CM  of  the  vehicle, 
o  inertial  acceleration  vector  of  the  CM  of  the  vehicle, 

o  attitude  of  the  vehicle  body  frame,  expressed  as  a  quaternion  and  dehned  as  a  3-2-1 
rotation  from  the  inertial  frame, 

o  angular  velocity  of  the  vehicle  body  frame  expressed  in  the  vehicle  body  frame, 
o  angular  acceleration  of  the  vehicle  body  frame  expressed  in  the  vehicle  body  frame,  and 
o  seconds  elapsed  since  the  beginning  of  the  trajectory. 

•  Initial  epoch,  which  corresponds  to  the  start  of  the  tabulated  truth  trajectory  and  is  specified 
as  a  Julian  Date. 
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•  Free  Navigation  Drift  time,  which  is  the  wait  time  for  the  vehicle  between  the  specified  initial 
epoch  and  the  beginning  of  the  trajectory  specified  by  the  tabulated  truth  data.  The  vehicle 
will  begin  navigating  at  the  initial  epoch,  and  is  assumed  fixed  to  the  Earth’s  surface  until 
the  Free  Navigation  Drift  time  has  elapsed. 

•  which  is  the  time  offset  between  Coordinated  Universal  Time  (UTC)  and  the  raw 
uncorrected  Universal  Time  (UT)  corrected  for  polar  wandering  (UTl). 

•  IMU  location  aboard  the  vehicle  with  respect  to  the  CM.  SAIMUN  synthesizes  the  centrifugal 
and  Euler  accelerations  experienced  by  the  IMU  due  to  the  angular  velocity  and  angular 
acceleration  of  the  vehicle. 

•  IMU  timestep,  which  is  the  time  between  successive  measurements  of  the  acceleration  and 
angular  velocity  provided  by  the  IMU. 

•  IMU  mechanization,  which  can  be  any  one  of  the  following: 

o  strap  down, 
o  space-stabilized,  or 
o  local-level. 

The  IMU  mechanizations  are  further  described  in  Section  3.4. 

•  Number  of  Monte  Carlo  sample  trajectories.  This  parameter  is  a  trade  between  computational 
effort  required  and  the  accuracy  of  the  results. 

•  Covariance  of  the  initial  position,  velocity,  and  attitude  of  the  vehicle,  corresponding  to  the 
initial  epoch. 

SAIMUN  then  implements  Monte  Carlo  analysis  to  propagate  the  initial  uncertainty  forward  in 
time  and  outputs  the  following  plots: 

•  position  standard  deviation  over  time, 

•  velocity  standard  deviation  over  time, 

•  attitude  standard  deviation  over  time, 

•  three-dimensional  terminal  miss  distance  histogram  of  the  Monte  Carlo  samples, 

•  two-dimensional  terminal  miss  distance  histogram  of  the  Monte  Carlo  samples  projected  in 
the  East-North  plane, 

•  propagation  error  of  the  trajectory  dead-reckoned  with  the  uncorrupted  acceleration  and  an¬ 
gular  velocity,  defined  as  the  propagated  truth,  with  respect  to  the  tabulated  true  trajectory 
input,  which  is  a  measure  of  the  numerical  degradation  due  to  discretization  of  the  continuous¬ 
time  equations  of  motion, 

•  three-dimensional  point  cloud  of  the  terminal  positions  of  the  Monte  Carlo  sample  trajectories 
in  the  East-North-Up  (ENU)  coordinate  frame, 

•  two-dimensional  point  cloud  of  the  terminal  positions  of  the  Monte  Carlo  sample  trajectories 
projected  in  the  East-North  plane. 
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•  two-dimensional  point  cloud  of  the  terminal  positions  of  the  Monte  Carlo  sample  trajectories 
projected  in  the  Up-East  plane, 

•  two-dimensional  point  cloud  of  the  terminal  positions  of  the  Monte  Carlo  sample  trajectories 
projected  in  the  North-Up  plane,  and 

•  three-dimensional  plot  of  the  true  and  propagated  truth  trajectories. 

SAIMUN  also  outputs  the  following  information: 

•  Circular  Error  Probable  (CEP),  which  is  defined  as  the  radius  of  a  circle  centered  about  the 
propagated  true  terminal  position  that  contains  half  of  the  East-North  projected  terminal 
positions  obtained  from  the  Monte  Carlo  samples, 

•  Spherical  Error  Probable  (SEP),  which  is  defined  as  the  radius  of  a  sphere  about  the  prop¬ 
agated  true  terminal  position  that  contains  half  of  the  terminal  positions  obtained  from  the 
Monte  Carlo  samples,  and 

•  Height  Error  Probable  (HEP),  which  is  defined  as  the  altitude  limits  (defining  an  interval) 
centered  about  the  propagated  true  terminal  position  that  contains  half  of  the  terminal  alti¬ 
tudes  obtained  from  the  Monte  Carlo  samples. 

3.2  IMU  Modeling 

The  acceleration  and  angular  velocity  measured  by  an  IMU  are  corrupted  by  a  variety  of  error 
sources.  SAIMUN  accounts  for  the  most  dominant  sources  of  error  by  specifying  a  probability 
density  function  (pdf)  for  the  error  parameters  to  define  each  error  source.  The  pdfs  implemented 
in  SAIMUN  are 

•  the  Dirac  distribution,  which  is  used  to  represent  a  constant  parameter  specified  solely  by  the 
mean, 

•  the  normal  distribution,  which  represents  a  Gaussian  distributed  parameter  that  is  solely 
defined  by  a  mean  and  variance  (or  standard  deviation),  and 

•  the  uniform  distribution,  which  represents  a  parameter  that  is  equiprobable  between  lower 
and  upper  bounds. 

These  distributions  are  graphically  represented  in  Fig.  (1). 

The  IMU  model  accounting  for  error  sources  is  given  for  the  accelerometers  and  gyroscopes  as 

0"m,k  =  +  -Ala  +  -Ma)  {I  Sa)  {cLk  +  f*a,0  +  +  'Wa,k)  ~  E  {ba,o}^  (1®) 

^m,k  =  <^Q,k  (^{1  +  -Alg  -|-  Mg)  (/  -|-  Sg)  {u)k  +  bgfl  +  t  g  ^k  T  §  ,k)  ~  E  )  (it) 

where 


ak  is  the  true  acceleration  experienced  by  the  IMU  in  the  IMU  frame  at  tk, 
baft  is  the  startup  bias  of  the  accelerometers, 

ba^k  is  the  bias  of  the  accelerometer  at  tk,  which  changes  due  to  bias  instability, 
Wa^k  is  the  thermo-mechanical  zero-mean  white  noise  present  in  the  accelerometers, 
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Figure  1.  Three  Distributions  Implemented  for  Error  Parameters  in  SAIMUN 

Sa  is  the  scale  factor  error  matrix  of  the  accelerometers, 

Ma  is  the  axes  misalignment  matrix  of  the  accelerometers, 

Na  is  the  axes  nonorthogonality  matrix  of  the  accelerometers,  and 

Uq  is  the  quantization  effect  caused  by  analog-to-digital  conversion. 

and  similarly  for  the  gyroscopes.  The  error  sources  will  now  be  presented  as  they  affect  the  ac¬ 
celerometer.  The  sources  of  error  used  by  this  model  affect  the  gyroscopes  identically  and  thus 
their  presentation  is  omitted. 

3.2.1  Startup  Bias.  The  startup  bias  affecting  the  accelerometers,  hafi-,  is  the  average  output  of 
the  accelerometers  while  they  are  not  undergoing  any  acceleration  [1].  For  practical  applications, 
the  startup  bias  is  quantified  by  accounting  for  and  removing  the  accelerometer  output  due  to 
gravity  while  the  IMU  remains  motionless.  Figure  (2)  shows  an  example  of  the  startup  bias  of  an 
accelerometer  undergoing  no  acceleration.  SAIMUN  assumes  that  the  statistics  of  the  startup  bias 


Time  [s] 


Figure  2.  Accelerometer  Output  Corrupted  by  Startup  Bias 

are  well  known  and  can  be  approximated  by  one  of  the  three  supported  error-parameter  pdfs.  The 
startup  bias  is  then  realized  from  the  error-parameter  pdf  during  virtual  IMU  initialization. 

3.2.2  Bias  Instability.  The  bias  of  an  accelerometer  tends  to  drift  slowly  over  time  and  is 
known  as  a  walking  bias  [1].  The  walking  bias  affecting  the  accelerometers  at  tk,  ba^k,  is  the  change 
in  the  bias  since  the  IMU  was  powered  on.  This  effect  is  modeled  by  a  first-order  random  walk 
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model,  which  causes  a  linear  increase  in  the  variance  of  the  accelerometer  bias  with  time.  The 
discrete  first-order  random  walk  model  is  given  by 


^a,k  —  ^a,k—l  T  ‘^Bl,k  ; 


(2) 


where  Wbi.A:  is  a  zero- mean  white  noise  process  of  standard  deviation  The  initial  standard 
deviation  in  the  walking  bias  is  zero,  and  the  standard  deviation  at  a  later  time,  tgj,  is  specified 
as  (Tgj.  Utilizing  the  zero  initial  condition  and  applying  Eq.  (2)  recursively  from  to  to  tk  yields  the 
walking  bias  at  as 


tklAt 

t^a,k  —  ^  ^  )  (3) 

i=l 

where  At  is  the  timestep  of  the  IMU.  Noting  that  the  variance  of  the  summation  of  n  independent 
random  variables  denoted  by  Xi  is  equivalent  to  the  sum  of  the  variances  of  n  independent  random 
variables,  i.e. 


this  yields 


Var 


E''”W ; 

i=l 


tk 

Var  [ba^k]  =  ^  Var  [wbi,*;]  =  ^<^1:  > 
i=l 


(4) 


when  applied  to  Eq.  (3).  Evaluating  Eq.  (4)  at  the  boundary  condition  gives 


tl  2  *  2 


(5) 


The  standard  deviation  of  the  zero- mean  white  noise  process  driving  the  random  walk  model,  <Tbi, 
can  now  be  found  by  manipulating  Eq.  (5)  as 


a 


BI  — 


(6) 


Eigure  (3)  shows  an  example  accelerometer  undergoing  a  random  walk  due  to  bias  instability. 
SAIMUN  inputs  the  boundary  conditions  <t*j  and  fgj,  the  IMU  timestep,  At,  and  the  desired 
error-parameter  pdf-type  for  and  calculates  the  standard  deviation  of  the  underlying  noise 
process  driving  the  random  walk  model  according  to  Eqs.  (2)-(6).  Each  step  in  the  random  walk 
model  is  realized  independently  for  each  IMU. 


3.2.3  Velocity  Random  Walk.  When  thermo-mechanical  zero-mean  white  noise,  fc,  present 
in  the  accelerometers  is  integrated,  a  random  walk  is  produced  in  the  velocity  state,  and  the 
standard  deviation  of  the  velocity  grows  proportionally  to  the  square-root  of  time  [1].  The  variance 
of  the  white  noise  present  in  the  accelerometers  is  described  by  the  Velocity  Random  Walk  (VRW) 
specification.  The  standard  deviation  of  the  zero-mean  white  noise  is  given  as  a  function  of  the 
VRW  specification  as 

_  VRWspec 

<^VRW  —  r-r - 

\/At 
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Figure  3.  Accelerometer  Output  Corrupted  by  Bias  Instability 


Time  [s] 


Figure  4.  Accelerometer  Output  Corrupted  by  Thermo-Mechanical  Zero-Mean  White 
Noise 


Example  accelerometer  output  corrupted  by  zero- mean  white-noise  is  shown  in  Fig.  (4).  SAIMUN  in¬ 
puts  the  VRW  specification  and  IMU  timestep  to  calculate  the  standard  deviation  of  the  white  noise 
process  underlying  the  velocity  random  walk.  This  zero-mean  white  noise  process  is  realized  inde¬ 
pendently  at  each  for  the  virtual  IMU. 


3.2.4  Scale  Factor  Error.  Errors  in  analog-to-digital  converters  and  other  circuitry  can  induce 
scale  factor  errors  in  the  accelerometers,  denoted  by  Sa-  For  instance,  errors  may  be  incurred  due 
to  reference  voltage  errors  caused  by  a  variety  of  sources,  including,  among  other  things,  tempera¬ 
ture  fluctuations  and  calibration  errors.  For  illustration  purposes,  consider  a  simple  example  of  a 
generic  scale  factor  error,  s,  and  assume  that  the  nominal  reference  voltage  for  the  analog-to-digital 
converter  is  Unom  =  5  V,  while  the  actual  reference  voltage  provided  to  the  analog-to-digital  con¬ 
verter  is  14ct  =  5.02  V.  The  scale  factor  error  of  the  signal  being  converted  by  the  analog-to-digital 
converter  is  then  given  by 


hiiom  U] 


S  = 


act 


act 


5  V  -  5.02  V 
5.02  V 


=  -0.003984  =  -0.3984%  . 


Denoting  the  scale  factor  errors  for  each  axis  of  the  accelerometer  triad  by  Sx,  Sy,  and  Sz,  the  total 
accelerometer  scale  factor  error  may  be  expressed  as 


Sa;  0  0 

0  Sy  0 

0  0  S2 


SAIMUN  inputs  an  error-parameter  pdf  for  Sx,  Sy,  Sz  and  realizes  these  parameters  for  the  virtual 
IMU  from  this  pdf. 
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3.2.5  Axes  Misalignment  Error.  Errors  in  the  manufacturing  process  and  mounting  of  an 
IMU  can  lead  to  axes  misalignment  errors,  denoted  by  Ma-  That  is,  the  actual  sensor  axes  are 
offset  from  the  IMU  case  frame  by  a  set  of  small  rotations.  This  is  readily  seen  by  considering  the 
term  I  +  Ma  in  Eq.  (la),  which  is  of  the  form  of  a  rotation  matrix  for  small  angles.  Therefore, 
the  actual  sensed  outputs  of  the  IMU  are  not  aligned  perfectly  with  the  modeled  case  frame  of  the 
IMU.  The  accelerometer  misalignment  errors  are  represented  within  Ma  by  nix,  niy,  and  m^,  with 
each  term  corresponding  to  its  respective  sensor  axis.  The  misalignment  errors  are  expressed  in 
matrix  form  as 


Ma 


0  ruz  —'niy 
—mz  0  nix 

niy  —nix  0 


SAIMUN  inputs  an  error-parameter  pdf  for  nix,  niy,  niz  and  realizes  these  parameters  for  the 
virtual  IMU  from  this  pdf. 


3.2.6  Axes  Nonorthogonality  Error.  Errors  in  the  manufacturing  process  of  an  IMU  can 
lead  to  axes  nonorthogonality  errors,  denoted  by  Na-  If  the  sensor  had  no  nonorthogonality  errors, 
each  axis  of  the  sensor  triad  would  be  at  perfect  right  angles  to  the  other  two  axes.  Therefore, 
the  nonorthogonality  errors  account  for  small  deviations  from  perfectly  orthogonal  axes  and  allow 
for  the  modeling  of  sensor  outputs  that  effectively  measure  small  overlaps  between  the  individual 
axes.  The  accelerometer  axes  nonorthogonality  errors  are  represented  within  Na  by  Ux,  ny,  and 
Uz,  with  each  term  corresponding  to  its  respective  sensor  axis.  The  nonorthogonality  errors  can  be 
expressed  in  matrix  form  as 


Na 


0 

Uz 

ny 

Uz 

0 

Ua 

.^y 

nx 

0 

SAIMUN  inputs  an  error-parameter  pdf  for  Ux,  Uy,  nz  and  realizes  these  parameters  for  the  virtual 
IMU  from  this  pdf. 

Example  output  of  an  accelerometer  triad  undergoing  1  m/s  acceleration  in  the  z  direction  with 
scale  factor,  axes  misalignment,  and  axes  nonorthogonality  errors  is  shown  in  Fig.  (5). 


3.2.7  Quantization  Error.  Quantization  error  arises  in  IMUs  due  to  the  limitations  of  analog- 
to-digital  signal  conversion  and  is  modeled  as  affecting  the  measurements  of  acceleration  and  angular 
velocity  as 

/2n-i 

aQ,k{ak)  =  round  (  - at 

V  ®max 

where  n  is  the  bitrate  of  the  analog  to  digital  converter,  and  Omax  is  the  maximum  value  able  to 
be  quantized  by  the  analog  to  digital  converter.  It  is  noted  that  the  quantization  error  is  applied 
after  all  of  the  other  errors  have  been  applied,  such  that  the  output  of  the  quantization  error  is 
the  measured  acceleration  of  Eq.  (la).  Values  of  the  bitrate  and  maximum  value  are  both  input 
into  SAIMUN.  The  quantization  error  in  angular  velocity  is  calculated  similarly.  Figure  (6)  shows 
a  quantized  signal  from  an  IMU  with  thermo- mechanical  noise  present. 


®max 
2n-l  ’ 
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Figure  5.  Three- Axis  Accelerometer  Output  Corrupted  by  Scale  Factor,  Misalignment, 
and  Nonorthogonality  Errors 


Time  [s] 


Fignre  6.  Accelerometer  Output  Corrupted  by  Quantization  Error 


3.3  Allan  Variance 

Section  3.2  outlined  the  modeling  of  the  dominant  error  sources  present  in  the  measured  acceleration 
(and  equivalently  the  angular  velocity)  signals  of  an  IMU.  In  order  to  utilize  these  models  the 
necessary  parameters  must  be  quantified  and  input.  One  technique  to  quantify  the  parameters 
necessary  to  model  the  velocity  random  walk  and  the  bias  instability  parameters  is  Allan  variance 
analysis  [1,  2,  3].  Allan  variance  analysis  quantifies  the  variance  between  averaging  times  of  a 
signal  in  order  to  determine  its  underlying  noise  process.  A  typical  Allan  variance  plot  is  shown  in 
Fig.  (7). 

The  slope  of  the  Allan  variance  curve  on  a  log-log  plot  corresponds  to  the  type  of  random 
process  dominating  the  signal  over  the  given  averaging  time.  The  slope,  m,  of  the  Allan  deviation 
(the  square  root  of  the  Allan  variance)  curve  for  different  underlying  noise  processes  present  in  an 
IMU  signal  consists  of  the  following; 

•  m  =  —0.5  corresponds  to  a  white  noise  process, 
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Figure  7.  Underlying  Noise  Processes  in  Allan  Variance  Analysis  [3] 


•  m  =  0  corresponds  to  bias  instability,  and 

•  m  =  0.5  corresponds  to  a  random  walk  process. 

To  illustrate  the  slopes  corresponding  to  each  of  the  aforementioned  noise  processes,  an  ex¬ 
ample  of  computing  the  Allan  deviation  for  synthetic  accelerometer  is  considered.  Three  different 
accelerometers  are  investigated:  a  baseline  accelerometer,  an  accelerometer  with  double  the  ve¬ 
locity  random  walk  specification  (relative  to  the  baseline),  and  an  accelerometer  with  double  the 
bias  instability  specification  (relative  to  the  baseline).  Data  on  each  accelerometer  is  simulated  for 
1  day  at  a  frequency  of  100  Hz.  This  data  is  then  used  to  construct  the  Allan  deviation  plot  given 
in  Fig.  (8)  using  the  algorithm  outlined  in  the  Appendix.  After  these  plots  are  generated  for  the 
accelerometers,  the  curves  can  be  used  to  extract  the  velocity  random  walk  and  bias  instability 
specifications  [1]. 

Comparison  of  the  baseline  and  double  VRW  Allan  deviation  curves  shows  that  an  increase 
in  white  noise  caused  by  the  increased  VRW  specification  causes  an  upward  vertical  shift  in  the 
portion  of  the  curves  representing  white  noise,  which  has  m  =  —0.5.  Because  the  bias  instability 
specifications  are  identical  for  the  baseline  and  double  VRW  curves,  the  portion  of  the  Allan 
deviation  curves  representing  the  random  walk,  which  has  m  =  0.5,  are  approximately  collinear. 

Comparison  of  the  baseline  and  double  bias  instability  Allan  deviation  curves  shows  that  an 
increase  in  bias  instability  causes  an  upward  vertical  shift  in  the  portion  of  the  curves  representing 
random  walk,  which  has  m  =  0.5.  Because  the  VRW  specifications  are  identical  for  the  baseline 
and  double  bias  instability  curves,  the  portion  of  the  Allan  deviation  curves  representing  the  white 
noise,  which  has  m  =  —0.5  are  approximately  collinear. 

The  velocity  random  walk  specification,  VRWgpec,  is  determined  by  the  intersection  of  a  line 
fitted  to  the  portion  of  the  Allan  deviation  curve  corresponding  to  the  white  noise  process  [m  = 
—0.5)  with  an  averaging  time  of  one  second.  The  intersection  point  determining  VRWgpec  for  the 
previously  described  baseline  accelerometer  is  illustrated  by  the  blue  marker  in  Fig.  (9).  Table  3 
summarizes  this  specification  for  the  baseline  accelerometer.  The  bias  instability  parameters, 
and  <T*j,  are  determined  by  the  minimum  point  of  the  Allan  deviation  curve.  This  point  is  shown  by 
the  green  marker  for  the  baseline  accelerometer  in  Fig.  (9),  and  the  specifications  are  summarized 
for  the  baseline  accelerometer  in  Table  1. 

Utilizing  Allan  variance  analysis  allows  for  the  measured  acceleration  and  angular  velocity  from 
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Allan  Deviation- Accelerometers  [m/s^]  Allan  Deviation- Accelerometers  [m/s^] 


Averaging  Time  [s] 

Figure  8.  Allan  Variance  Analysis  of  Example  Accelerometers 


Averaging  Time  [s] 

Figure  9.  Allan  Variance  Analysis  for  the  Baseline  Accelerometer 
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Table  1.  VRW  and  Bias  Instability  Specifications  from  Allan  Variance  Analysis  for 
the  Baseline  Accelerometer 


Specification 

Value 

VRWspec 

9.694  X  10-4 

26.56  s 

2.679  X  10-4  ^ 

different  IMUs  to  be  compared  and  contrasted  through  trends  in  the  Allan  deviation  curves.  While 
the  above  discussion  has  focused  on  Allan  variance  analysis  as  applied  to  accelerometers,  the  same 
process  can  be  used  when  considering  gyroscopes. 

3.4  IMU  Mechanizations 

Traditionally,  IMU-based  navigation  is  mechanized  in  one  of  two  ways:  gimbaled/gyro-stabilized  or 
strapdown.  In  the  gimbaled  mechanization,  the  triad  of  accelerometers  is  placed  on  a  “stable  table” 
at  the  center  of  set  of  interconnected  gimbals.  This  enables  the  stable  table,  through  the  use  of 
feedback  tracking  using  the  gyros,  to  track  a  predetermined  attitude  such  that  the  accelerometers 
measure  a  certain,  desirable  resultant  linear  motion.  Typically,  the  stable  table  is  configured  to 
track  a  north-east-down  coordinate  system,  in  which  case  the  accelerometers  measure  the  linear 
motion  in  these  directions,  or  it  is  configured  to  track  the  inertial  coordinate  system  (i.e.  to  remain 
fixed  in  space  as  the  vehicle  maneuvers). 

The  most  basic  of  gimbal-based  mechanizations  employs  three  interconnected  gimbal  rings,  such 
as  is  found  on  the  Minuteman  III  intercontinental  ballistic  missile  (see  Fig.  (10(a))).  While  this 
configuration  provides  the  ability  to  actively  control  the  orientation  of  the  stable  table  through 
motors  mounted  on  the  gimbals,  the  use  of  three  gimbals  sometimes  leads  to  gimbal  lock  and 
prevents  the  stable  table  from  tracking  all  possible  attitudes  that  may  be  required. 


(a)  Minuteman  III  IMU 


(b)  Peacekeeper  IMU 


Figure  10.  Examples  of  Stable  Table  IMU  Mechanizations  [4] 

More  advanced  conhgurations  of  the  stable  table  mounting  design  can  be  used  to  replace  the 
original  gimbaled  system.  This  is  accomplished  by  using  four  or  more  gimbals  or  by  using  a  fluid 
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stabilization  system,  which  has  the  stable  table  floating  in  the  center  of  a  sphere  filled  with  pressur¬ 
ized  fluid.  Whereas  the  attitude  of  the  stable  table  is  controlled  using  torquers  and  manipulation 
of  the  gimbals  in  the  original  design,  the  attitude  of  the  stable  table  in  the  fluidic  design  is  con¬ 
trolled  via  nozzle  actuators  which  avoids  any  possible  complications  due  to  gimbal  lock.  One  such 
implementation  of  the  fluidic  system  is  the  Advanced  Inertial  Reference  Sphere  (AIRS)  that  was 
developed  for  the  LGM-118A  Peacekeeper  (see  Fig.  (10(b))). 

In  the  strapdown  mechanization,  the  triads  of  accelerometers  and  gyros  are  fixed  (strapped 
down)  to  the  body  of  the  vehicle,  and  there  is  no  stable  table.  This  means  that  where  the  stable 
table  was  used  to  mechanically  track  the  inertial  frame  in  order  to  track  a  predetermined  attitude 
and  sense  accelerations  in  this  rotating  frame,  the  gyro  outputs  are  now  used  to  numerically  track 
the  orientation  of  the  IMU  (called  the  case  frame)  with  respect  to  the  inertial  frame.  To  avoid  any 
singularities  present  in  three-parameter  attitude  representations  such  as  Euler  angles,  the  attitude 
quaternion  is  used  and  can  be  propagated  forward  in  time  using  the  sensed  angular  velocity  from 
the  IMU.  Additionally,  the  accelerations  are  now  sensed  in  the  case  frame,  and  these  must  be 
transformed  using  the  case-to-inertial  attitude  in  order  to  predict  the  translational  position  and 
velocity  of  the  vehicle. 

Three  inertial  navigation  system  mechanizations  are  presented  here;  a  strapdown  system,  an 
actively  commanded  local-level  gimbaled  system,  and  a  passively  commanded  inertially  stable  space- 
stabilized  gimbaled  system.  For  each  of  the  aforementioned  mechanizations,  the  governing  equations 
for  the  position,  velocity,  and  attitude  of  the  IMU  are  presented,  and  a  set  of  discretized  equations 
are  developed.  The  discretized  equations  governing  the  position,  velocity,  and  attitude  of  the  IMU 
then  are  representative  of  the  software  and  hardware  that  would  be  used  in  each  inertial  navigation 
mechanization. 

3.4.1  Strapdown  Mechanization.  The  inertial  position,  r,  velocity,  v,  and  attitude,  q,  of  a 
strapdown  IMU  aboard  a  vehicle  are  governed  by  [5,  6] 

=  «g(uVu(^)  +  ^c(^)^CM/lMu)  + 

^i{t)  = 

where  a  is  acceleration,  ca  is  angular  velocity,  and  is  a  rotation  matrix  that  rotates  a  vector 
representation  from  the  c  to  the  i  coordinate  frame.  Q  represents  the  angular  velocity  expressed  as 
a  pure  quaternion.  The  c/i  subscript  notation  represents  the  c  frame  with  respect  to  the  i  frame. 
The  sub/superscripts  refer  to 

IMU  inertial  measurement  unit 
CM  center-of-mass  of  the  vehicle 
c  case  frame  of  the  IMU 
i  Earth  Centered  Inertial  frame 
g  due  to  gravity 
ng  not  due  to  gravity 

where  the  superscript  for  a  vector  quantity  denotes  the  coordinate  frame  in  which  the  vector  is 
expressed.  To  simplify  the  nomenclature,  let 


UMu(i) 

r{t)  , 

XX)  - 

v{t)  , 

X) 

9{-) 

TX 

T^{t)  , 

P 

1 

a-{t)  > 

Qtit) 

Qit) 

u:{t)  , 

/mC  _ 

cm/imu 

-)■  d , 

and  rl^X)  +  TX)X/,Mv 

s{t) 
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With  these  substitutions,  the  equations  of  motion  may  be  written  succinctly  as 


r{t)  =  v{t) 

v{t)  =  g{s{t))  +T'^{t)a{t) 

The  continuous  time  strapdown  IMU  equations  of  motion  are  now  discretized  via  analytical 
integration.  In  the  development  of  these  discretized  equations  of  motion,  it  is  assumed  that  the 
non-gravitational  acceleration  and  angular  velocity  are  constant  over  a  small  time-step,  such  that 

f^m,k  —  vT  and  i^m,k  —  yv  • 

IXtk 

The  continuous  time  attitude  evolves  as 

^{t)  =  (8)  q{t) , 

where  u)m(t)  is  the  pure  quaternion  given  (in  right-handed,  vector- first  form)  by 


Um{t) 


0 


Define  the  quaternion  Aq{t)  representing  the  rotation  from  the  a  priori  attitude  as 


Aq{t)  =  q{t)  (g)  . 


Then,  solving  for  q{t)  yields 

q{t)  =  Aq{t)  (g)  Qk-i  ■ 

Assuming  the  time  step  to  be  small,  it  can  be  shown  that  the  rotation  vector  from  tk-i  to  tk  is 

^{tk)  —  ^m,k  {tk  tk—l)  =  ^m^kAtk  , 

where  Atk  =  tk  —  tfc-i-  Thus,  the  discrete  quaternion  propagation  is  given  by 

qk  =  q{tk)  =  q{^m,kAtk)  ^  Qk-i ,  (7) 

where  the  quaternion  representation  of  the  arbitrary  rotation  vector  0  is  given  by 

sin  (^6*)  0/9 


m  = 


cos 


The  continuous  time  velocity  evolves  as 

v{t)  =  g{s{t))  +  T^{t)am{t) . 


Assume  a  small  time  step  such  that  the  non-gravitational  acceleration  is  approximately  constant, 
i.e.  ttmit)  =  (im,k-  Then,  it  follows  that 

v{t)  =  g{s{t)) +  T'^{t)a^^k-  (8) 
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The  gravitational  acceleration  may  be  expanded  in  a  first-order  Taylor  series  expansion  about  the 
a  priori,  yielding 

g{s{t))  =  g{sk-i)  +  Gk-i  {s{t)  -  Sk-i)  , 

where 


^fc-i  —  f‘k-1  + 


and 


Gk-i  = 


dg{s) 

ds 

From  the  relationship  that  s  =  r  +  T^d,  and  letting  r{t)  — )■  Vk-i,  it  follows  that 

g{s{t))  =  gisk-i)  +  Gk-i  d. 

The  first-order  evolution  of  the  rotation  matrix  is  given  by 

T{t)  =  Tk-i  —  [6{t)x]Tk-i , 


(9) 


which  may  be  transposed  to  give 


T^(t)=Tj_i+Tj_i[0(t)x].  (10) 

Applying  this  relationship  to  the  gravitational  acceleration  in  Eq.  (9),  noting  that  6{t)  =  — 

tk-i),  and  defining  gk-i  =  g{sk-i),  yields 


g{s{t))  =  gk-i  -  Gk-iT^_i[dx]u:)rn,k  {t  -  4-i)  •  (11) 

Substituting  Eqs.  (10)  and  (11)  into  Eq.  (8),  simplifying,  and  rearranging  terms  yields 

v{t)  =  gk-i  +  T^_iam,k  -  {Gk-iT^-i[dx]  [am,k^])  i^m,k  {t  -  4-i)  •  (12) 

Integrating  Eq.  (12)  twice  and  once,  respectively,  yields  the  position  and  velocity  evolution  in  time 
as 

r{t)  =  Tk-i  +  Vk-i  {t  -  tk-i)  +  ^gk-i  {t  -  tk-if  +  ]^T^_iam,k  (t  -  tk-if 

~  g  {Gk-lT^_i[dx]  +  T^-i  \o,m,k  x])  (^  ~  ^fc-l)^ 

v{t)  =  Vk-I  +  gk-i  {t  -  tk-i)  +  T^_iam.,k  {t  -  tk-i) 

~  2  {^k-lTk-l[dx]  +  Tj[_i  [am,k  x])  <^m,k  {t  —  tk-l)"^  ■ 

Letting  t  — )■  tfc  yields  the  discrete  propagation  equations  for  the  position  and  velocity  of  the  IMU 
as 


Tk  =  Vk-i  +  t’fc-iAffc  -F  ]^T'^_^am,k^tl 


(13a) 


~  g  {Gk-lTk-i[dx\  -t-  T'^_i  [am,fcX])  LJm,k^tk 

Vk  =  Vk-i  +gk-iAtk  +  T^_iam,kAtk  -  ^  {Gk-iT^_i[dx]  +  [a^.fcx])  Um,kAtl .  (13b) 

In  summary,  the  discrete-time  dead-reckoning  equations  for  the  position,  velocity,  and  attitude 
of  a  vehicle  with  a  strapdown  IMU  are  given  by  Eqs.  (13)  and  (7). 
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3.4.2  Local-Level  Mechanization.  The  local-level  mechanization  physically  commands  the 
IMU  platform  via  torquers  into  alignment  with  the  local  geographic  frame  tracking  north,  east, 
and  down  (NED).  This  class  of  navigation  system  maintains  two  of  its  accelerometers  along  the 
local  horizon  with  the  third  aligned  normal  to  the  reference  ellipsoid.  To  keep  this  alignment,  the 
system  rotates  at  an  angular  velocity  equal  to  the  Earth’s  rotation  rate  and  the  vehicle  motion 
along  the  Earth  and  is  independent  of  the  vehicle  body  attitude.  The  navigational  quantities  of 
interest  are  the  IMU  frame  position  and  Earth-referenced  velocity  vectors.  A  representation  of  the 
NED  coordinate  frame  is  shown  in  Fig.  (11). 

^ecef 


Figure  11.  Representation  of  the  NED  Coordinate  Frame^ 


The  geodetic  latitude,  longitude,  and  negative  altitude  are  governed  by  [7]: 


1 


A  = 


{ri  +  h) 

1 


In  -  {tl  +  ‘^h)^  -  ^{ri  +  h){X  +  2uJie)Xsm{2(f>)  +  ig 


Ie 

{ri  -|-  /i)  cos  4> 


+  {tl  +  U  +  2h)(l)Xiaii(f)  +  2{rL  +  h)4)Uiet&'a.(j) 


-  2h[X  +  L0ie)  -  h  'x  - 


Jk_ 

cos  (j) 


-h  =  Id  -  in  +  h){X  +  2uJie)Xcos^  4>  -  {ri  +  +  g, 


(14a) 

(14b) 


(14c) 


where  cj)  is  the  geodetic  latitude,  A  is  the  longitude,  h  is  the  altitude,  Wje  is  the  scalar  rotation  rate 
of  the  Earth  about  its  spin  axis.  The  terms  /v,  /e,  and  fn  are  components  of  the  specific  force 
vector,  /sj,  which  describe  the  values  of  the  specific  force  in  the  north,  east,  and  down  directions 
respectively.  These  are  the  values  of  the  sensed  IMU  acceleration  in  the  local-level  mechanization. 

The  magnitude  of  the  geocentric  position  vector  to  the  point  on  the  ellipsoid  directly  below  the 
vehicle  is  denoted  by  tq.  The  meridional  radius  of  curvature,  r^,  is  the  ellipsoid  radius  of  curvature 

^  Image  modified  with  permission  from 

http:/ /en.wikipedia.org/wiki/File:ECEF_ENU_Longitude_Latitude-relationships.svg 
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in  the  north-south  direction  and  the  prime  radius  of  curvature,  ri,  is  normal  to  the  meridional 
radius  of  curvature  in  the  east-west  direction  on  the  ellipsoid.  The  equations  for  computing  ro,  vl, 
and  ri  are 

ro  =  re  [l  -  ^/(l  -  cos(2(/>))  -h  -  cos(4(/)))] 

=  'To  [l  —  2/ cos(2(/))  -|-  |/^(1  —  4 cos(2(/))  -|-  7 cos(4(/)))] 
r;  =  ro  [l  -|-  2/  sin^  (j)  —  3/^  sin^  (/>(!  —  2  sin^  i^)]  . 

From  the  equations  for  rg,  rr,  and  ri,  temporal  derivatives  of  the  meridional  radius  of  curvature 
and  the  prime  radius  of  curvature  can  be  determined.  Furthermore,  the  evaluation  of  the  radii 
depends  upon  Ve  and  /,  which  are,  respectively,  the  equatorial  radius  and  the  flattening  of  the 
Earth  ellipsoid.  The  flattening  factor,  /,  is  a  measure  of  the  compression  between  the  equatorial 
and  polar  radii  and  is  defined  to  be 


where  is  the  polar  radius  of  the  Earth.  The  three  quantities  rg,  r^,  and  /  are  constant  for 
the  Earth  reference  ellipsoid;  for  instance,  when  using  the  1984  World  Geodetic  System  (WGS-84) 
model  for  the  reference  ellipsoid,  the  parameters  are  given  as  shown  in  Table  2  (note  that  Vp  is  a 
derived  constant). 

Table  2.  Parameters  Defining  the  WGS-84  Model  (adapted  from  [8]) 


Parameter 

Value 

re 

6378.1370  km 

rp 

6356.7523  km 

/ 

1/298.257223563 

^ie 

7.2921150  X  10-5  rad/s 

3.986004418  x  10^  km^/s^ 

The  final  elements  required  to  complete  the  evaluation  of  Eqs.  (14)  are  the  components  of  the 
gravity  vector,  ^g,  rjg,  and  g.  These  are  the  components  of  the  gravity  vector  referenced  in  the 
geographic  (local-level)  frame,  which  is  given  by 


1 

_ 1 

4 

Gtv  —  rojf^  sin  4>  cos  4>c 

9^  = 

1 

1 _ 

— 

-Vg 

_ge  +  A5(_ 

— 

Ge 

Gd  —  ruf^  cos  4>  cos  4>c 

where  (/>c  is  the  geocentric  latitude,  ^g  is  the  gravitational  component  due  to  meridian  deflection 
of  the  vertical,  gg  is  the  gravitation  acceleration  due  to  prime  deflection  of  the  vertical,  is  the 
gravity  magnitude  associated  with  the  reference  ellipsoid,  and  A.g  is  the  gravity  anomaly,  which  is 
the  difference  between  the  observed  gravitational  acceleration  and  the  value  predicted  from  a  given 
model.  Values  of  the  gravity  anomaly  are  available  online  from  the  National  Geospatial-Intelligence 
Agency.  The  geocentric  latitude,  i;f)c,  is  defined  as  the  angle  between  the  equatorial  plane  and  the 
normal  and  is  given  by  the  difference  between  the  geodetic  latitude  and  the  deviation  of  the  normal 


</>c  =  </>  -  /(I  -  h/ro)  sin  {24)) . 
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To  complete  the  evaluation  of  ^g,  r]g,  and  g,  Gn,  Ge,  and  Gd  are  approximated  as 

Gn  ~  -G^  -  Grfsm{2(f)) 

Ge  =  0 

Ge  ~  —Gr  +  Gfj,f  sin(2(/>) , 

where  Gr  and  are  the  radial  and  colatitude  components  of  the  gravitational  field  in  spherical 
coordinates.  Retaining  oblateness  effects  up  to  and  including  the  fourth  order  terms,  Gr  and  G^f, 
may  be  written  as 


where  ip  is  the  vehicle’s  colatitude. 

The  discrete-time  form  of  the  latitude,  longitude,  and  negative  altitude  dynamics  is  formed  using 
an  Euler-integration  method  on  the  continuous-time  dynamics.  Defining  a  “position,”  “velocity,” 
and  “acceleration”  to  be 

</>  1  [01 

r  =  A  ,  V  =  X  ,  and  a  =  A  , 

_—h_  _—h_  _—h_ 

the  discrete  propagation  equations  are  given  by 

Tk  =  Tk-i  +  Vk-lAtk 

Vk  —  Vk—\  -|-  f^k—lAtk  , 

where  (^,  A,  and  —h  are  given  by  the  evaluation  of  Eqs.  (14). 

In  the  local-level  IMU  mechanization,  the  stable  table  tracks  the  local-level  (NED)  coordinate 
system,  which  changes  orientation  as  the  vehicle’s  latitude  and  longitude  change.  The  orientation 
of  the  local-level  coordinate  frame  in  reference  to  the  inertial  frame  is  given  by  the  direction  cosine 
matrix 

—  sin  (j)  cos  A 
T^=  —  sin  (j)  sin  A 

cos  4> 

where  A  is  the  celestial  longitude,  given  by 

A  =  A  —  Aq  -|-  (jJiet  , 

where  A  is  the  terrestrial  longitude  from  Greenwich,  Aq  is  the  initial  terrestrial  longitude,  and  t  is 
the  elapsed  time.  The  angular  velocity  of  the  local- level  frame  is 

A  cos  4> 

-XsinKp 


—  sin  A  —  cos  (j)  cos  A 
cos  A  —  cos  4>  sin  A  , 
0  —  sin  (j) 


(15a) 

(15b) 
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which  may  be  used  to  form  the  time-wise  evolution  of  the  orientation  of  the  local-level  frame  via  [5] 

^hit)  =  ®  Quit)  , 

where  is  the  pure  quaternion  (in  right-handed,  vector-first)  form  of  the  local-level  frame 

angular  velocity,  defined  as 


For  notational  ease,  let  —S'  and  qh{t)  — )■  q{t). 

The  discrete  evolution  of  the  IMU  platform  attitude  is  as  follows.  Define  the  quaternion  Aq{t) 
representing  the  rotation  from  the  a  priori  attitude  as 

Aqit)  =  q{t)  <S)  q^li  ■ 


Solving  for  q{t)  yields 


q{t)  =  Aq{t)  (g)  Qfc-i  • 

Assuming  the  time  step  to  be  small,  it  can  be  shown  that  the  rotation  vector  from  tk-i  to  tk  is 


6{tk)  =  ujk  {tk  —  tk-i)  =  ^kAtk  , 

where  Atk  =  tk  —  tfc-i  and  =  u)^^^{tk)-  Thus,  the  discrete  quaternion  propagation  is  given  by 


qk  =  q{tk)  =  q{^kAtk)  ®  qk-i , 


(16) 


where  the  quaternion  representation  of  the  arbitrary  rotation  vector  0  is  given  by 

sin  (lO)  0/9 
cos (le) 


In  summary,  the  discrete  time  evolution  of  the  local-level  position,  velocity,  and  IMU  platform 
attitude  are  governed  by  Eqs.  (15)  and  (16). 


3.4.3  Space-Stabilized  Mechanization.  The  space-stabilized  mechanization  tracks  an  Earth- 
centered,  inertially  nonrotating  coordinate  frame.  This  reference  frame  is  realized  by  a  gyro- 
stabilized  platform  with  triaxial  control.  Except  for  torques  applied  as  compensation  to  account 
for  anisoelastic  effects,  thermal  and  acceleration  sensitivity,  etc.,  the  gyroscopes  are  uncommanded 
and  the  system  is  “free-floating.” 

The  Earth-centered  inertial  acceleration  of  the  IMU  is  taken  as  [7] 

r\t)  =  fjt)+g\r^it)), 

where  r*  is  the  inertial  position  of  the  IMU,  is  the  inertial  gravitational  acceleration  experienced 
by  the  IMU,  and  accelerometer  misalignment  and  nonorthogonality  errors  have  been  absorbed  into 
the  measured  specific  force  vector  /^.  It  should  be  noted  that  any  displacement  effects  caused 
by  the  gravitational  acceleration  acting  at  the  CM  as  opposed  to  the  IMU  are  neglected.  Noting 
that  the  specific  force  is  actually  measured  in  the  IMU  platform  frame,  the  translational  equations 
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governing  the  dynamical  evolution  of  the  position  of  the  IMU  can  be  expressed  as  a  set  of  first-order 
differential  equations  given  by 


r*(t)  =  v^{t)  (17a) 

e{t)  =  g\r\t))  +  Ti,{t)fm,  (17b) 

where  is  the  inertial  velocity  of  the  IMU,  fm  is  the  sensed  acceleration  in  the  IMU  platform 
frame,  and  is  the  attitude  transformation  from  the  platform  frame  to  the  inertial  frame.  Ideally, 
the  platform  frame  exactly  tracks  the  inertial  frame;  if  this  is  the  case,  then  Tp  =  /sxs-  The 
continuous-time  attitude  representation  of  the  platform  frame  evolves  as  [5] 

^iit)  =  ,  (18) 

where  q{t)  is  the  pure  quaternion  (in  right-handed,  vector-first)  form  of  the  IMU  frame  angular 
velocity,  defined  as 


pit 


it) 


For  notational  simplicity,  let  r*(f)  — ^  r(t),  v^{t)  — ^  v{t),  ^^p/^it)  — >  <^mit),  qfit)  — >  qit),  Tp{t)  = 
and  fmit)  =  fm{t),  such  that  the  equations  governing  the  forward  evolution  of  the  position, 
velocity,  and  attitude  of  the  space-stabilized  IMU  are  given  by  collecting  Eqs.  (17)  and  (18)  and 
making  the  appropriate  substitutions  to  arrive  at 


r{t)  =  v{t) 

vit)  =girit))  +  T^{t)fmit) 

Ht)  =  ^^mit)  (8)  qit) . 

The  continuous-time  equations  of  motion  are  now  discretized  via  analytical  integration.  In  the 
following  equations,  it  is  assumed  that  the  non-gravitational  acceleration  and  angular  velocity  are 
constant  over  the  time-step: 

1  ^^m,k 

Tm,k  —  W7  ana  ^rn,k  —  TT  • 

The  continuous  time  attitude  evolves  as 

Ht)  =  ^^mit)  (8)  qit) , 

where  i^mit)  is  the  pure  quaternion  given  (in  right-handed,  vector- first  form)  by 


LJmit) 


LJmit) 

0 


Define  the  quaternion  Aq{t)  representing  the  rotation  from  the  a  priori  attitude  as 

Aqit)  =  qit)  (8)  q^l^ . 


Solving  for  q{t)  yields 


q{t)  =  Aq{t)  (8)  Qfc-i  • 
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Assuming  the  time  step  to  be  small,  it  can  be  shown  that  the  rotation  vector  from  to  is 


^{tk)  —  ^m,k{tk  tk—l)  —  ^m,k^tk  ; 

where  At^  =  —  tfc-i-  Thus,  the  discrete  quaternion  propagation  is  given  by 

Qk  =  q{tk)  =  ®  Qfc-i , 

where  the  quaternion  representation  of  the  arbitrary  rotation  vector  6  is  given  by 


(19) 


q{e)  = 


sin  (^0)  G/d 
cos (i0) 


The  continuous  time  velocity  evolves  as 

v{t)  =  g{r{t))  +  T^{t)fm{t) 


Assuming  a  small  time  step  such  that  the  non- gravitational  acceleration  is  approximately  constant, 
i.e.  am{t)  =  such  that 

v{t)  =  gir{t))  +  T'^{t)fm,k  ■ 

The  first-order  evolution  of  the  rotation  matrix  is  given  by 

T{t)  =T{q{t))  =  Tfc_i  -  [e{t)x]Tk.i , 
which  may  be  inverted  to  give 


T^(t)=Tj_i+Tj_i[0(t)x]. 

Letting  r{t)  Vk-i  and  defining  g{r{t))  =  g{rk-i)  gk-i  yields 

V{t)  =  {Tk_ifm,k  +  fl'fc-i)  {t  -  tk-l)  . 

Integrating  twice  and  once  respectively,  and  letting  t  ^  tk  gives  the  following; 

Tk  =  Tk-i  +  Vk-iAtk  +  ^  {Tk-ifm,k  +  gk-i)  (20a) 

Vk  =  Vk-i  +  {T^_ifm,k  +  gk-i)  ^tk  ■  (20b) 

In  summary,  the  inertial  position,  velocity,  and  IMU  platform  attitude  of  a  space-stabilized 
inertial  navigation  system  evolve  in  discrete  time  according  to  Eqs.  (20)  and  (19). 

3.5  Monte  Carlo  Analysis 

In  general,  any  of  the  mechanizations  of  the  IMU  can  be  represented  by  a  nonlinear  dynamical 
system.  If  the  relevant  expressions  for  position  (e.g.  latitude,  longitude  and  altitude  or  inertial 
Cartesian  position),  velocity,  and  attitude  are  collected  into  a  state  vector,  x,  the  general  nonlinear 
dynamical  system  representing  the  forward  evolution  of  the  state  for  any  of  the  IMU  systems  may 
be  expressed  as 


—  f  (®fc— 1 )  tk)  j 


(21) 
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where  and  cj^  represent  the  true  vehicle  non-gravitational  acceleration  and  angular  velocity,  re¬ 
spectively.  For  the  mechanizations  presented  in  Section  3.4,  /(•)  represents  one  of  Eqs.  (13)  and  (7), 
Eqs.  (15)  and  (16),  or  Eqs.  (20)  and  (19).  Eunctionally,  the  relationship  between  the  true  and  mea¬ 
sured  IMU  outputs  is  given  by 

am,k  =  a{4>k,tk)  and  u^^k  =  u{il)k,tk) ,  (22) 

where  (pk  and  ipk  represent  the  IMU  model  parameters,  which  also  have  governing  evolutionary 
equations  as 

(pk  =  g{(pk,  tk)  and  ipk  =  Hipk,  tk)  ■ 

In  the  preceding  discussions,  many  possibilities  for  the  forms  of  g,  h,  cp,  and  ip  have  been  given.  For 
example,  for  the  simple  additive  model  where  the  true  IMU  outputs  are  corrupted  by  a  constant 
bias  b  and  a  white- noise  sequence  w,  Eqs.  (22)  become 

^m,k  —  4“  ‘^a,k  and  ^m,k  —  ^k  4“  bg  -|-  Wg^k  ; 

and  the  parameters  of  the  IMU  error  model  are 


4>k  — 

ba 

and 

-ipk  = 

_  Wa,k  _ 

.  _ 

with  the  governing  equations  cpk  =  (pk-i  and  ipk  =  ipk-i  (note  that  the  white  noise  sequence  does 
not  require  a  formal  time- wise  evolution  owing  to  the  properties  of  white  noise). 

In  order  to  predict  the  position,  velocity,  and  attitude  of  the  vehicle,  the  standard  approach  is 
to  implement  an  approximate  filter  such  as  the  extended  Kalman  filter  (EKF).  The  EKE  is  a  linear 
minimum  mean  square  error  filter  which  does  not  require  an  assumption  of  Gaussianity  except  to 
make  statistical  conclusions.  In  the  EKF,  the  state  estimate  is  represented  by  x,  and  the  IMU 
outputs  are  used  to  formulate  a  state  estimate  evolution  of  the  form 

—  f  {^k—l:  tk)  ■> 

which  is  accompanied  by  a  linearized  propagation  of  the  uncertainty  distribution  of  the  states 
through  the  covariance  propagation  equation 

Pk  =  Fk-iPk-iFk-i  +  Qk  ■  (23) 

In  Eq.  (23),  Fk-i  is  the  tangent  linear  dynamics  evaluated  at  the  current  estimated  state  and 
estimated  IMU  model  parameters,  and  Qfc  is  the  process  noise  covariance  matrix.  It  should  be 
noted  that  the  linearization  in  the  EKE  is  performed  about  the  current  best  state  estimate.  The  two 
preceding  equations  then  represent  the  best  estimate  of  the  vehicle  state  along  with  the  uncertainty 
in  that  estimate. 

Typical  analysis  of  IMU-driven  systems  employs  linearization  about  a  reference  truth  trajectory 
in  order  to  decouple  the  covariance  propagation  from  the  state  estimate  propagation;  this  is  what  is 
called  linear  covariance  analysis.  In  doing  so,  the  true  IMU  acceleration  and  angular  velocity  that 
are  used  in  the  forward  solution  to  Eq.  (21)  are  also  used  in  the  covariance  propagation  equation  of 
Eq.  (23),  obfuscating  the  appearance  of  any  influence  from  the  IMU  errors  on  the  state  estimate  and 
covariance.  Furthermore,  the  covariance  obtained  from  the  linear  covariance  analysis  is  interpreted 
as  being  representative  of  a  Gaussian  distribution  in  order  to  determine  performance  measures; 
however,  if  the  distribution  is  not  Gaussian,  the  calculation  of  performance  measures  can  lead 
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to  erroneous  conclusions  when  using  the  assumption  of  Gaussianity.  Additionally,  the  mean  and 
covariance  information  from  a  Gaussian  distribution  gives  no  realistic  representation  of  the  closest 
approach  of  an  estimate  to  the  truth  or  the  furthest  departure  of  an  estimate  from  the  truth;  that 
is,  no  best  case  or  worst  case  scenario  information  can  be  obtained.  In  order  to  obviate  erroneous 
statistical  conclusions  based  on  a  Gaussian  assumption  resulting  from  using  the  covariance  matrix 
obtained  from  a  linear  covariance  analysis,  Monte  Carlo  techniques  that  statistically  sample  the 
IMU  error  parameterization,  simulate  sensed  IMU  outputs,  emulate  the  navigation  filter  in  order 
to  obtain  estimates  of  the  vehicle’s  position,  velocity,  and  attitude,  and  then  collect  statistics  on 
certain  user-defined  are  needed. 

To  illustrate  the  Monte  Carlo  simulation  process,  consider  the  two  dimensional  transformation 
from  polar  to  Cartesian  coordinates,  i.e. 


X 

y 


Given  a  distribution  on  m  =  [r  0],  it  is  desired  to  determine  distribution  characteristics  on 
z  =  [x  y].  First,  assume  that  r  and  9  are  independent  and  Gaussian  with  mean  and  covariance 


1  [ml 

r  0.02  [ml^  0  1 

rriw  = 

L  J 

60  [deg] 

and  Pw  = 

L  J 

0  30  [deg]^ 

Therefore,  the  polar  to  Cartesian  conversion  may  be  written  as 


2:  =  f{w) , 

where  w  is  the  input  and  z  is  the  associated  output.  The  Monte  Carlo  simulation  process  is 
started  by  drawing  a  number  of  samples  from  the  input  vector,  where  the  sampling  is  performed 
with  respect  to  the  probability  density  function  which  governs  the  characteristics  of  the  input. 
Once  the  samples  are  drawn,  the  nonlinear  mapping  function,  f,  is  applied  to  each  input  sample 
to  generate  a  set  of  output  samples. 

For  the  coordinate  conversion  problem,  the  Monte  Carlo  simulation  is  initialized  by  drawing 
1  X  10^  realizations  of  the  random  vector  w  from  the  distribution 

=  |27rPi„|"^/^exp  |-^  {w  -  m^)'^  P~^  {w  -  m^) 

The  set  of  samples  drawn  is  shown  in  Fig.  (12(a)).  The  nonlinear  transformation  is  applied  to  each 
sample  in  the  input  space  and  the  resulting  set  of  Cartesian  coordinates  is  shown  in  Fig.  (12(b)).  It 
is  clear  that  the  output  space  is  not  capable  of  being  described  by  a  Gaussian  distribution;  however, 
the  Monte  Carlo  simulation  technique  can  still  be  used  to  compute  an  accurate  approximation  of 
the  mean  and  covariance  via 

1  ^ 

and  P2  =  —  ^(zj -m2)(zj -m^)'^,  (24) 

i=l 


1 


N 


i=l 


where 


z*  =  f{Wi) 

and  Wi  is  the  Monte  Carlo  sample.  The  contours  representing  the  Gaussian  distribution  com¬ 
puted  from  the  preceding  mean  and  covariance  equations  are  shown  on  top  of  the  set  of  output-space 
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(a)  Input  Space  Monte  Carlo  Samples  (b)  Mapped  Samples  in  the  Output  Space 

Figure  12.  Input  and  Ontput  Samples  from  the  Monte  Carlo  Simulation 


Monte  Carlo  samples  in  Fig.  (13(a)).  Moreover,  a  linearization  approach  can  be  used  to  map  the 
input  mean  and  covariance  into  an  approximate  output  mean  and  covariance.  The  comparison 
between  the  linear  mean  and  covariance  and  the  Monte  Carlo  mean  and  covariance  is  shown  in 
Fig.  (13(b)),  wherein  it  is  seen  that  the  two  differ  a  great  deal.  The  only  characteristic  that  is 
mostly  captured  by  the  linearization  process  is  the  orientation  of  the  covariance  matrix.  The  cen¬ 
tering  (i.e.  the  mean  of  the  distribution)  is  not  captured  accurately,  and  the  size  of  the  mapped 
covariance  matrix  is  not  captured  accurately. 


(a)  Monte  Carlo  Mean  and  Covariance 


(b)  Comparison  between  Linearization  and  Monte  Carlo 


Figure  13.  Distribution  Contours  Obtained  by  the  Mean  and  Covariance  from  the 
Monte  Carlo  Simulation  and  Linearization 


To  further  demonstrate  the  power  of  Monte  Carlo  simulation,  consider  a  modification  of  the 
preceding  problem.  In  this  case,  it  is  still  assumed  that  the  radial  and  angular  input  coordinates 
are  independent,  but  now  the  radial  direction  is  taken  to  be  Gaussian  distributed  and  the  angular 
direction  is  taken  to  be  uniformly  distributed.  The  mean  and  covariance  are  taken  to  be  the  same 
as  before  so  that  the  input  space  has  the  same  first  and  second  central  moments  as  in  the  preceding 
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example.  The  probability  density  function  for  the  input  is  now  given  by 


p{w,mr,Pr,m0,Pe) 


1 271  Pr  I  exp 


{Wr  —  nirY  \  1 

J  ’  2^/Wb 

{Wr  —  rrir)'^  '[ 

J  ’ 


where  mo  —  \/3Pg  <  6  <  mo  +  \/3Po,  such  that  [mo  —  \/3Po  ,  mo  +  \/3Po]  is  the  support  set  of  the 
uniform  angular  distribution.  From  the  original  specifications  for  the  transformation  problem,  the 
means  and  covariances  required  to  evaluate  the  input-space  distribution  are 


mr  =  I  [m]  ,  me  =  60  [deg] ,  Pr  =  0.02  [m]^  ,  and  Po  =  30  [deg]^  . 

As  before,  the  Monte  Carlo  simulation  is  initialized  by  drawing  1  x  10^  realizations  of  the  random 
vector  w  from  the  distribution  p{w;mr,  Pr,mg,  Po),  and  each  sample  is  mapped  into  the  output 
space  via  the  nonlinear  transformation.  These  samples  are  shown  in  Figs.  (14(a))  and  (14(b)), 
respectively.  The  mean  and  covariance  is  again  computed  from  the  Monte  Carlo  samples  and 
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(b)  Mapped  Samples  in  the  Output  Space 


Figure  14.  Input  and  Output  Samples  from  the  Monte  Carlo  Simulation  (modified 
problem) 

shown  along  with  the  output-space  samples  in  Fig.  (15(a)).  Finally,  a  linearization  approach  is 
used  to  map  the  input-space  mean  and  covariance  into  an  output-space  mean  and  covariance. 
The  comparison  between  the  Monte  Carlo  mean  and  covariance  and  the  linearization  mean  and 
covariance  is  shown  in  Fig.  (15(b)). 

While  it  is  clear  that  the  linearization  process  once  again  fails  to  capture  the  true  mean  and 
covariance,  it  is  noted  that  the  linearized  mean  and  covariance  depicted  in  Fig.  (15(b))  are  identical 
to  the  ones  depicted  in  Fig.  (13(b)).  At  the  same  time,  however,  the  Monte  Carlo  means  and 
covariances  differ  between  these  two  simulations.  This  occurs  because  the  Monte  Carlo  simulation 
has  and  uses  knowledge  of  the  entire  distribution  on  both  the  input  and  output  spaces  whereas  a 
linearization  procedure  works  with  only  the  first  and  second  central  moments  of  the  distribution. 
Therefore  if  two  input  distributions  have  the  same  mean  and  covariance,  and  are  subjected  to  the 
same  nonlinear  transformation,  the  linearized  solutions  will  always  have  the  same  output  mean  and 
covariance. 
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1.2 


(a)  Monte  Carlo  mean  and  covariance 


(b)  Comparison  between  linearization  and  Monte  Carlo 


Figure  15.  Distribution  Contours  Obtained  by  the  Mean  and  Covariance  from  the 
Monte  Carlo  Simnlation  and  Linearization  (modified  problem) 


It  is  also  worth  noting  that  the  Monte  Carlo  simulation  enables  analysis  extending  beyond  the 
computation  of  the  mean  and  covariance.  For  instance,  higher-order  statistical  moments  may  be 
computed  in  a  similar  fashion  to  the  mean  and  covariance,  best  case  and  worst  case  performance 
can  be  analyzed  through  analysis  of  the  dispersions  away  from  a  predefined  truth,  or  the  radial 
distance  encapsulating  a  given  percentage  of  the  sample  points  can  be  computed. 

The  problem  of  analyzing  the  performance  of  IMU-based  navigation  systems  is  characterized 
by  complex,  nonlinear  functions  where  the  input  space  is  described  by  the  IMU  error  sources, 
which  may  have  fairly  arbitrary  statistical  descriptions.  Monte  Carlo  analysis  naturally  handles 
the  characteristics  of  this  problem  and  enables  a  more  rich  data  analysis  suite  than  is  found  in 
typical  covariance-based  analysis. 

4.0  RESULTS  AND  DISCUSSION 

4.1  Trajectory  Simulation 

In  order  to  provide  the  trajectory  input  to  SAIMUN  that  is  outlined  in  Section  3.1,  a  ground  ve¬ 
hicle  was  used  to  log  acceleration  and  angular  velocity  data  using  a  VectorNav  VN-100^  IMU.  The 
VN-100  logged  data  at  40  Hz  during  a  trajectory  that  took  approximately  41  minutes  to  complete. 
A  smoother  was  implemented  to  decrease  the  effect  of  noise  on  the  recorded  acceleration  and  angu¬ 
lar  velocity  data.  The  smoothed  acceleration  and  angular  velocity  data  were  dead-reckoned  from 
a  known  starting  condition  to  provide  the  true  trajectory  for  SAIMUN  to  use.  The  dead-reckoned 
trajectory  is  not  perfectly  representative  of  the  trajectory  followed  by  the  vehicle;  however,  for 
the  purposes  of  applying  SAIMUN  and  analyzing  IMU-based  navigation  performance,  the  dead- 
reckoned  trajectory  is  taken  to  be  the  true  trajectory. 

4.2  IMU  Selection 

To  demonstrate  the  utility  of  SAIMUN,  two  inertial  measurement  units  are  compared:  the  Lord 
Microstrain  3DM-GX3-15^  and  the  Analog  Devices  ADIS16488^.  Both  IMUs  utilize  MEMS  sensors 

^http:/ /www.vectornav.com/Downloads/Support/PB-12-0002.pdf,  accessed  04/21/2014 
^http:/ /files. microstram.com/3DM-GX3-15-Inertial-Measurement-Unit. pdf,  accessed  04/21/2014 
"^http:/ /www. analog.com/static/imported-files/data_sheets/ADIS16488. pdf,  accessed  04/21/2014 
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to  provide  triaxial  accelerometers,  gyroscopes,  magnetometers,  and  a  pressure  sensor.  The  IMUs 
are  implemented  to  provide  dead- reckoning  navigation  of  a  ground  vehicle  trajectory. 

For  this  analysis,  only  the  accelerometer  and  gyroscopic  information  is  of  interest.  Table  3  shows 
the  pertinent  specifications  for  both  of  the  chosen  IMUs.  These  IMUs  are  interesting  to  compare 
due  to  the  differences  between  the  random  walk  and  bias  instability  specifications.  Table  3  shows 
that  the  ADIS  16488  has  a  higher  random  walk  specification  and  thus  more  white  noise  present  in 
the  signal.  The  3DM-GX3-15,  on  the  other  hand,  has  a  higher  bias  instability  specification.  The 
higher  bias  instability  and  lower  random  walk  specifications  of  the  3DM-GX3-15  compared  to  the 
ADIS16488  imply  that  the  position  and  attitude  uncertainty  associated  with  the  3DM-GX3-15  will 
grow  more  slowly  than  the  ADIS16488  initially.  Eventually,  the  position  and  attitude  uncertainty 
associated  with  the  3DM-GX3-15  will  grow  larger  and  surpass  the  position  and  attitude  uncertainty 
associated  with  the  ADIS16488  as  time  increases. 

Table  3.  IMU  Specifications  for  the  ADIS16488  and  3DM-GX3-15 


Specification 

3DM-GX3-15 

ADIS16488 

Analog-to-Digital  Converter  Bitrate 

16 

32 

Range 

±156.96  f 

±176.58  f 

C/2 

o; 

Scale  Factor  Error 

^f(-500,500)  ppm 

7f(-5000,5000)  ppm 

o; 

a 

Axes  Nonorthogonality  Error 

^(—103,103)  arcsec 

7/(— 126, 126)  arcsec 

o 

Axes  Misalignment  Error 

^(—103,103)  arcsec 

7/(— 3600,  3600)  arcsec 

o 

o 

Velocity  Random  Walk 

0.00078481^ 

vnr 

0.029  1^ 

vnr 

C 

Bias  Instability 

0.0003924  f 

0.000981 

Bias  Instability  Time 

100  s 

100  s 

Startup  Bias 

(-0.0196, 0.0196)  f 

7f(-0.016,  0.016)  f 

Analog-to-Digital  Converter  Bitrate 

16 

32 

Range 

±600  f 

±450  f 

cc 

Scale  Factor  Error 

7f(-500,500)  ppm 

7f(-100,100)  ppm 

0) 

Cl. 

O 

Axes  Nonorthogonality  Error 

7/(— 103, 103)  arcsec 

7/(— 180, 180)  arcsec 

O 

cc 

O 

Axes  Misalignment  Error 

7/(— 103, 103)  arcsec 

7/(— 3600,  3600)  arcsec 

o 

Angular  Random  Walk 

0.03  ^ 

vnr 

Bias  Instability 

18  s 

6.25^ 

Bias  Instability  Time 

100  s 

100  s 

Startup  Bias 

7f(-0.25,0.25)  f 

U{-0.2,0.2)  1 

4.3  Performance  Metrics 

In  order  to  assess  the  performance  of  an  IMU  for  a  given  trajectory,  metrics  need  to  be  defined. 
These  metrics  include  error  probables,  miss  distance  histograms,  and  standard  deviation  in  the 
position,  velocity,  and  attitude.  Because  the  ADIS16488  and  the  3DM-GX3-15  are  both  relatively 
low-end  IMUs,  two  cases  are  considered:  Case  1  analyzes  only  the  first  500  seconds  of  the  entire 
41  minute  trajectory  while  Case  2  analyzes  the  entire  41  minute  trajectory.  For  each  performance 
metric  considered,  the  analysis  is  carried  out  for  each  of  the  two  cases  that  are  considered. 
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4.3.1  Error  Probables.  The  CEP  represents  the  radius  of  a  circle  centered  on  the  true  terminal 
position  in  the  East-North  plane  in  which  the  probability  of  landing  in  the  circle  is  0.5.  The  CEP 
is  typically  based  on  a  circular  bivariate  normal  distribution,  which  relies  on  the  assumption  that 
the  terminal  points  of  all  samples  will  follow  a  Gaussian  distribution.  SAIMUN,  however,  does  not 
use  the  typical  Gaussian  assumption  to  calculate  the  CEP.  Since  Monte  Carlo  sampling  techniques 
are  implemented,  the  CEP  is  calculated  by  finding  the  radius  under  which  50%  of  the  East-North 
plane  projection  of  terminal  points  of  the  Monte  Carlo  samples  fall. 

The  SEP  represents  the  radius  of  a  sphere  centered  on  the  true  terminal  position  in  which  the 
probability  of  landing  in  the  sphere  is  0.5.  Again,  since  Monte  Carlo  sampling  is  used,  no  Gaussian 
assumption  is  used  to  calculate  the  SEP. 

The  HEP  represents  the  altitude  range  centered  on  the  true  terminal  position  in  which  the 
probability  of  landing  within  this  range,  whether  it  is  above  or  below  the  ground,  is  0.5.  Again,  no 
Gaussian  assumption  is  used  to  calculate  the  HEP  due  to  the  advantage  of  Monte  Carlo  sampling. 

All  error  probables  use  the  true  trajectory  at  the  hnal  time  to  calculate  the  terminal  position. 
The  error  probables  for  the  ADIS16488  and  the  3DM-GX3-15  are  shown  in  Table  4  for  both  Case  1 
and  Case  2.  The  values  of  the  error  probable  are  computed  by  analyzing  the  terminal  Monte  Carlo 
samples.  For  example,  when  considering  the  SEP,  the  three-dimensional  miss  distance  is  computed 
for  each  sample.  These  values  are  then  sorted,  and  the  radius  which  encompasses  50%  of  the  miss 
distances,  which  is  the  SEP,  is  determined. 

Table  4.  Error  Probable  Values  for  the  ADIS16488  and  3DM-GX3-15 


Error  Probable 

MicroStrain  3DM-GX3-15 

Analog  Devices  ADIS  16488 

t-H 

CEP 

1.527  km 

1.674  km 

§ 

SEP 

1.536  km 

1.719  km 

u 

HEP 

0.132  km 

0.319  km 

CN 

CEP 

2757  km 

1000  km 

§ 

SEP 

2862  km 

1114  km 

u 

HEP 

665  km 

364  km 

4.3.2  Miss  Distance  and  Terminal  Error.  The  miss  distance  represents  the  magnitude  of 
the  terminal  position  error  with  respect  to  the  propagated  true  trajectory  at  the  final  time.  Two 
metrics  are  considered;  a  two-dimensional  miss  distance  and  a  three-dimensional  miss  distance. 
The  two-dimensional  metric  uses  the  magnitude  of  the  East  and  North  terminal  position  error 
components,  and  the  three-dimensional  miss  distance  is  calculated  using  the  magnitude  of  the 
full  ENU  frame  error  vector.  By  calculating  the  miss  distance  of  each  sample  trajectory  that  is 
generated  by  the  Monte  Carlo  simulation  and  using  a  histogram  plot,  the  sample  error  distribution 
is  determined.  The  two-dimensional  and  three-dimensional  miss  distance  histograms  are  given  in 
Fig.  (17(a))  for  the  Case  1  analysis  and  in  Fig.  (17(b))  for  the  Case  2  analysis. 

ENU  projections  of  the  terminal  position  error  for  every  Monte  Carlo  sample  are  shown  in 
Figs.  (18)”(21).  Figure  (18)  illustrates  the  projections  for  the  ADIS16488  IMU  as  applied  to  the 
Case  1  trajectory.  Fig.  (19)  illustrates  the  projections  for  the  3DM-GX3-15  IMU  as  applied  to 
the  Case  1  trajectory.  Fig.  (20)  illustrates  the  projections  for  the  ADIS16488  IMU  as  applied 
to  the  Case  2  trajectory,  and  Fig.  (21)  illustrates  the  projections  for  the  3DM-GX3-15  IMU  as 
applied  to  the  Case  2  trajectory.  It  is  worth  noting  that  the  axis  limits  for  the  same  projections 
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Figure  16.  Miss  Distance  Histograms  for  Two-  and  Three-Dimensional  Positions:  Case 
1 


in  Figs.  (18)  and  (19)  are  identical,  such  that  a  direct  comparison  on  the  size  of  the  uncertainty 
can  be  made.  Similarly,  this  is  true  for  Figs.  (20)  and  (21).  Figures  (18)-(21)  illustrate  the 
non-Gaussian  nature  of  the  underlying  pdf.  Conventional  uncertainty  analysis  techniques,  such 
as  linear  covariance  analysis,  which  only  capture  the  first  and  second  statistical  moments  of  the 
underlying  pdf,  rely  on  an  assumption  to  be  made  about  the  underlying  pdf  in  order  to  calculate 
performance  metrics  such  as  CEP,  SEP,  and  HEP.  This  assumption  of  the  underlying  pdf,  usually 
that  it  is  Gaussian,  can  lead  to  unrepresentative  performance  metrics  to  be  calculated.  Monte  Carlo 
analysis  allows  the  assumption  of  an  underlying  pdf  to  be  relaxed  because  the  mean,  covariance, 
CEP,  SEP,  and  HEP  can  be  calculated  from  the  discrete  Monte  Carlo  samples. 

4.3.3  Standard  Deviation  in  Position  and  Attitude.  In  order  to  quantify  uncertainty  in  the 
position  and  attitude  states,  the  standard  deviation,  or  square  root  of  the  second-central  moment, 
is  calculated  and  plotted  as  a  function  time  based  on  the  Monte  Carlo  samples  (see  Eqs.  (24)). 
While  beyond  the  scope  of  the  current  analysis,  the  standard  deviation  as  computed  from  the 
Monte  Carlo  samples  can  be  compared  to  other  methods  of  uncertainty  propagation,  such  as  linear 
covariance  analysis.  The  position  and  attitude  standard  deviations  computed  for  the  ADIS16488 
and  3DM-GX3-15  IMUs  can  be  found  in  Eig.  (24)  for  the  Case  1  trajectory  and  in  Eig.  (25)  for  the 
Case  2  trajectory. 
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Figure  17.  Miss  Distance  Histograms  for  Two-  and  Three-Dimensional  Positions:  Case 
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Figure  18.  Terminal  Position  Error  of  the  Monte  Carlo  Samples  for  the  ADIS16488 
IMU:  Case  1 
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Figure  19.  Terminal  Position  Error  of  the  Monte  Carlo  Samples  for  the  3DM-GX3-15 
IMU:  Case  12 
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Figure  20.  Terminal  Position  Error  of  the  Monte  Carlo  Samples  for  the  ADIS16488 
IMU:  Case  2 
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Figure  21.  Terminal  Position  Error  of  the  Monte  Carlo  Samples  for  the  3DM-GX3-15 
IMU:  Case  2 
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Figure  22.  Position  Standard  Deviations  as  a  Function  of  Time  for  the  ADIS16488 
and  the  3DM-GX3-15  IMUs:  Case  1 
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Figure  23.  Attitude  Standard  Deviations  as  a  Function  of  Time  for  the  ADIS16488 
and  the  3DM-GX3-15  IMUs:  Case  1 
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Figure  24.  Position  Standard  Deviations  as  a  Function  of  Time  for  the  ADIS  16488 
and  the  3DM-GX3-15  IMUs:  Case  2 
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Figure  25.  Attitude  Standard  Deviations  as  a  Function  of  Time  for  the  ADIS  16488 
and  the  3DM-GX3-15  IMUs:  Case  2 
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5.0  CONCLUSIONS 


Before  choosing  an  IMU  to  be  implemented  in  a  navigation  system,  SAIMUN  can  be  used  to 
analyze  the  performance  of  the  IMU  along  a  given  trajectory  or  trajectories.  When  comparing 
the  Lord  MicroStrain  3DM-GX3-15  and  the  Analog  Devices  ADIS16488  IMUs  in  the  navigation 
of  an  example  vehicle  trajectory,  the  lower  random  walk  specification  of  the  3DM-GX3-15  causes 
its  associated  position  uncertainty  to  grow  slower  during  the  beginning  of  the  trajectory;  however, 
the  higher  bias  instability  of  the  3DM-GX3-15  causes  its  associated  position  uncertainty  to  grow 
faster  and  overtake  the  position  uncertainty  associated  with  the  ADIS16488.  For  these  IMUs, 
480-500  seconds  is  observed  to  be  the  approximate  window  in  which  the  performance  of  the  IMUs 
is  similar,  as  quantified  by  the  position  variance.  Before  this  window,  the  3DM-GX3-15  provides 
a  more  accurate  navigation  solution.  After  this  window,  the  ADIS16488  provides  a  more  accurate 
navigation  solution.  SAIMUN  allows  this  trade  study  to  be  performed  before  purchasing  or  testing 
any  hardware  and  will  allow  for  larger  trade  spaces  to  be  investigated  virtually,  obtaining,  building, 
or  testing  any  physical  IMUs. 

SAIMUN  enables  both  gimbal-  and  strapdown-based  navigation  systems  to  be  analyzed  virtually 
in  the  developed  simulation  environment.  An  IMU  error  parameterization  and  a  trajectory  are 
defined  in  the  simulation  environment  and  Monte  Carlo  sample  trajectories  are  dead- reckoned 
according  to  the  appropriate  dynamics  for  the  user-dehned  mechanization  type.  SAIMUN  negates 
the  need  for  any  assumption  to  be  made  about  the  underlying  distribution  of  the  navigation  solution 
to  calculate  performance  metrics,  such  as  is  necessary  in  the  more  traditional  linear  covariance 
analysis.  The  Monte  Carlo  analysis  implemented  in  SAIMUN  can  more  accurately  quantify  desired 
performance  metrics,  mean,  and  covariance,  without  requiring  any  underlying  assumptions  on  the 
distribution  of  the  state  uncertainties  to  be  made. 

Currently,  SAIMUN  is  capable  of  rendering  navigation  solutions  for  gimbal-  and  strapdown- 
based  navigation  systems,  while  considering  IMU  error  sources  due  to  startup  bias,  bias  instabil¬ 
ity,  thermo-mechanical  zero-mean  white  noise,  scale  factor  errors,  axes  misalignment  errors,  axes 
nonorthogonality  errors,  and  quantization  effects  caused  by  analog-to-digital  conversion.  Where 
appropriate,  SAIMUN  allows  for  a  selection  from  multiple  probabilistic  representations  of  the  IMU 
error  sources,  such  as  Dirac,  Gaussian,  and  uniform  distributions.  The  structure  of  SAIMUN  sup¬ 
ports  future  extensions,  such  that  additional  IMU  mechanizations,  IMU  parameter  errors,  and  error 
source  distributions  can  be  incorporated.  On  the  analysis  side  of  SAIMUN,  this  is  also  true;  that  is, 
analyses  extending  beyond  those  presented  in  this  report  can  be  incorporated.  The  simulation  and 
analysis  environment  is  designed  to  be  modular  such  that  these  modifications  are  straightforward 
to  implement  within  SAIMUN. 
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APPENDIX  -  CALCULATING  ALLAN  VARIANCE 


The  Allan  deviation  plot  of  generic  process  which  is  sampled  at  a  constant  frequency  fs,  can 
be  produced  by  the  following  algorithm; 

1.  Log  a  number  of  samples,  N,  from  the  generic  process  u}{t) 

N  =  fsT, 


where  T  is  duration  of  the  data  series. 

2.  Create  a  vector  of  monotonically  increasing  averaging  times,  defined  as 

taw  —  [f-l  U  •••  ^n]  ) 

where  ti  >  1/ fs  and  tn  <  T/2. 

3.  Divide  uj{t)  into  Mi  clusters  for  the  taw  as 

T 

Mi=  -  , 

ltaw,i\ 

where  f  =  1,  2,  ...,  n  and  the  [-J  function  represents  rounding  towards  negative  infinity. 

4.  For  each  cluster  Mj,  compute  the  time  average  of  Lo{t)  for  each  portion  of  cluster  size  Mi 


j=fsti,v,i{k-l)+l 


where  k  =  1,  2,  ...,  Mi 

5.  The  Allan  variance  for  the  averaging  time,  cr‘^{taw,i)  is  then  calculated  as 

^  M,-l 

O'  {taw,i)  =  2(^M-  —  1)  ~  tOi^k)  ■ 


6.  The  Allan  deviation,  a{taw,i),  is  then  plotted  as  a  function  of  the  averaging  time,  taw,i,  to 
produce  the  Allan  deviation  plot. 
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LIST  OF  SYMBOLS,  ABBREVIATIONS,  AND  ACRONYMS 


AIRS 

CEP 

CM 

EKE 

ENU 

HEP 

IMU 

MEMS 

NED 

pdf 

SAIMUN 

SEP 

UT 

UTl 

UTC 

VN 

VRW 

WGS-84 


Advanced  Inertial  Reference  Sphere 

Circular  Error  Probable 

Center  of  Mass 

Extended  Kalman  Eilter 

East-North-Up 

Height  Error  Probable 

Inertial  Measurement  Unit 

MicroElectricalMechanical  Systems 

North-East-Down 

Probability  Density  Punction 

Simulation  and  Analysis  of  IMU-based  Navigation 
Spherical  Error  Probable 
Universal  Time 

Universal  Time  corrected  for  polar  wandering 

Coordinated  Universal  Time 

VectorNav 

Velocity  Random  Walk 

1984  World  Geodetic  System 
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